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We have extended the density matrix renormalization group (DMRG) approach to two-fluid open 
many-fermion systems governed by complex-symmetric Hamiltonians. The applications are carried 
out for three- and four-nucleon (proton-neutron) systems within the Gamow Shell Model (GSM) 
in the complex-energy plane. We study necessary and sufficient conditions for the GSM-I-DMRG 
method to yield the correct ground state eigenvalue and discuss different truncation schemes within 
DMRG. The proposed approach will enable configuration interaction studies of weakly-bound and 
unbound strongly interacting complex systems which, because of a prohibitively large size of Fock 
space, cannot be treated by means of the direct diagonalization. 
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I. INTRODUCTION 



The theoretical description of weakly bound and un- 
bound states in atomic nuclei requires a rigorous treat- 
ment of many-body correlations in the presence of scat- 
tering continuum and decay channels P, [li Many- 
body states located close to the particle emission thresh- 
old display unusual properties such as, e.g., halo and Bor- 
romean structures, clusterization phenomena, and cusps 
in various observables resulting from the strong coupling 
to the continuum space. These peculiar features cannot 
be described in the standard shell model (SM) in which 
the single-particle (s.p.) basis is usually derived from 
an infinite well such as the harmonic oscillator poten- 
tial. Indeed, resonance and scattering states, which play 
a decisive role in the structure of weakly bound/unbound 
states, are not properly treated within the standard SM 
formalism. 

The solution of the configuration-interaction problem 
in the presence of continuum states has been recentl 
advanced in the real-energy continuum SM 0, [^, 
and in the complex-energy SM, the so-called Gamow 
SheU Model (GSM) 0, i, i, [HI- In the GSM, the 
s.p. basis is given by the Berggren ensemble containing 
Gamow states and the non-resonant continuum of scat- 
tering states. The scattering states are distributed along 
a contour defined in the complex fc-plane and, together 
with the Gamow states, form a complete set [l^]- In prac- 
tice, the contour is discrctized and the many-body basis 
is spanned by the Berggren ensemble. As in standard 
SM applications, the dimension of the many-body valence 
space increases dramatically with the number of valence 



nucleons and the size of the s.p. basis. Moreover, the use 
of the Berggren ensemble implies complex-symmetric ma- 
trices for the representation of the Hermitian Hamilton 
operator. Consequently, efhcient numerical methods are 
needed to solve the many-body Schrodinger equation of 
GSM. The DMRG approach is ideally suited to optimize 
the size of the scattering space in the GSM problem as 
the properties of the non-resonant shells vary smoothly 
along the scattering contour. 

The DMRG method was first introduced to over- 
come the limitations of the Wilson-type renormaliza- 
tion group to describe strongly correlated ID lattice sys- 
tems with short-range interactions [l^ (see recent re- 
views [3, m, [ill)- More recently, by reformulating 
the DMRG in a s.p. basis, several applications to fi- 
nite Fermi systems like molecules [l7j_ 18 1, superconduct- 
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clei [2^, and fractional quantum Hall systems [24J| have 
been reported. While most of the DMRG studies were 
focused on equilibrium properties in strongly correlated 
closed quantum systems characterized by Hermitian den- 
sity matrices, non-equilibrium systems involving non- 
Hermitian and non-symmetric density matrices can also 
be treated [1^. Nuclear applications of DMRG in the 
context of the standard SM, both in the M-scheme and in 
the angular-momentum conserving J-scheme, have also 
been reported [26.,^., 28| with mixed success. In the pre- 
vious study we reported the first apphcation of the 
DMRG method in the context of the GSM and showed 
that in this case the method provides a highly accurate 
description of broad resonances in the neutron-rich nuclei 
with few valence particles. 
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The present study is an extension of the previous work 
[29| to the case where both protons and neutrons are in- 
cluded in the valence space. Several significant improve- 
ments over Ref. have been made concerning the 
DMRG algorithm and numerical implementation. Our 
work is organized as follows. In Sec. HIl we briefly reca- 
pitulate the GSM formalism and the generalized vari- 
ational principle behind it. Section IIIII describes the 
DMRG method in the J-scheme as applied to the open- 
system formalism of the GSM. Practical applications of 
the GSM+DMRG method are presented in Sec.|TVl Illus- 
trative calculations are carried out for ''Li (three-nucleon 
systems) and *Li (four-nucleon systems). We study the 
impact of different starting conditions on the DMRG re- 
sult and introduce two different truncation schemes in the 
DMRG procedure and discuss their virtues. The trunca- 
tion schemes are compared with the GSM benchmark di- 
agonalization results. The resulting efficient calculation 
scheme opens a window for extending calculations to sys- 
tems beyond the current limits of direct diagonalization. 
The conclusions of our work are contained in Sec. fVl 
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FIG. 1: Illustration of the Berggren completeness relation ([T]) 
in the complex fc-plane. The bound states are located on the 
positive imaginary axis. The weakly bound halo states lie 
close to the origin. The positive-energy resonant states are 
located in the fourth quadrant. Those with a small imagi- 
nary part can be interpreted as resonances. The complex-fc 
shells on the contour L'^ represent the non-resonant scattering 
continuum. 



II. THE GAMOW SHELL MODEL 

The s.p. basis used in the GSM formalism is generated 
by a finite depth potential. It forms a complete set of 
states in the sense of the Berggren completeness relation 
El: 



\Un){Un\ + / 



\uk){uk\ = 1, 



(1) 
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where the discrete sum includes s.p. bound (negative- 
energy) resonant states (6) and positive-energy decaying 
resonant states (d) . These states are the poles of the cor- 
responding one-body scattering matrix. The integration 
in Eq. ([T]) is performed along a contour defined in the 
complex fc-plane that is located below the resonant states 
included in the basis (see Fig. [T|). In general, different 
contours can be used for each £, j partial wave. By dis- 
cretizing continuum states on i+ a finite s.p. basis is 
obtained. The many-body basis is obtained in the usual 
way by constructing product states (Slater determinants) 
from this discrete s.p. set. We assume in the follow- 
ing that the nucleus can be described as a system of 
protons and neutrons evolving around a closed core. 
Within this picture, the GSM Hamiltonian H reads: 



H = 



E 

i=l 
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where p?/2/i is the s.p. kinetic energy operator, /i is 
the reduced mass of the nucleon-|-core system, Ui is the 
finite-depth, one-body potential, and Vij is the two-body 
residual interaction. 

The GSM Hamiltonian ^ is complex symmetric. Ac- 
cording to the complex variational principle [sol . [sij . a 



complex analog to the usual variational principle for self- 
conjugated Hamiltonians, the Rayleigh quotient 



is stationary around any eigenstate |<I>o) of H: 
H\%) = E[%]\%). 



(3) 



(4) 



That is, at |$) = |$o) the variation of the functional 
is zero: 



= 0. 



(5) 



It should be noted that the complex variational 

principle is a stationary principle rather than an upper or 
lower bound for either the real or imaginary part of the 
complex eigenvalue. However, it can be very useful when 
applied to the squared modulus of the complex eigenvalue 
Indeed, 

S^lEl"^ = S^{E*E) = E*S^E + ES^E* ^0 (6) 
at |<f>) = |<f>o) because of analyticity of E[^]. 



III. THE DENSITY MATRIX 
RENORMALIZATION GROUP METHOD FOR 
THE GAMOW SHELL MODEL 

Let us consider a nucleus with active protons and rii, 
active neutrons and let us denote by \J^) the eigenstate 
of H having angular momentum J and parity tt. As \J^) 
is the many-body pole of the scattering matrix of H, the 
contribution from scattering shells on L"*" to the many- 
body wave function is usually smaller than the contribu- 
tion from the resonant orbits. Based on this observation, 
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the following separation is usually performed [29|: the 
many-body states constructed from the s.p. poles form 
a subspace A (the so-called 'reference subspace'), and 
the remaining states containing contributions from non- 
resonant shells form a complement subspace B. As we 
shall discuss later in Sections IIV Al and IIVB| this intu- 
itive definition of the reference subspace may be insufh- 
cient for describing certain classes of eigenstates. In such 
cases, states in A have to be constructed both from the 
s.p. poles and selected scattering shells. 

At the first stage of the GSM-f DMRG method, called 
'the warm-up phase', the scattering shells are gradually 
added to the reference subspace to create the subspace 
B. This process is described in the next section. 

A. Warm-up phase of GSM+DMRG 

One begins by constructing all states \k)A forming the 
reference subspace A. The set of those states shall be 
denoted as {kA}- The many-body configurations in A 
can be classified in different families {^i; j^} according to 
their number of nucleons n, total angular momentum j^, 
and parity tt. In the following, we shall omit the parity 
label in the notation of a given family. States with a 
number of protons (neutrons) larger than n,r (n,^) are 
not considered since they do not contribute to the many- 
body states in the composition of subspaces A and B. 

All possible matrix elements of suboperators of the 
two-body Hamiltonian ([2]) acting in A, expressed in the 
second quantization form, are calculated and stored: 

{O} = {a\{a^af,{a^a^f,{{a^a^)''a)^, 

{a^aYiaaf}, (7) 

with and a being the nucleon creation and annihilation 
operators in resonant shells. The GSM Hamiltonian is 
then diagonalized in the reference space to provide the 
zeroth-order approximation l^'j)'-"-' to \ J'^). This vector, 
called 'reference state', plays an important role in the 
GSM-I-DMRG truncation algorithm. 

In the next step, the subspace of the first scatter- 
ing shell {lj)i belonging to the discretized contour L+ 
is added. Within this shell, one constructs all possible 
many-body states {(0)i^}: denoted as grouped in 
{"-b; js} families. Matrix elements of suboperators ([7]) 
acting on \i)B are computed. By coupling states in A 
with the states lis), one constructs the set {/c^® is}"^ of 
states having fixed J'^. This ensemble serves as a basis in 
which the GSM Hamiltonian is diagonalized. Of course 
at this stage, the resulting wave function is a rather poor 
approximation of \J'^) as only one scattering shell has 
been included. The target state is selected among 
the eigenstates of H as the one having the largest overlap 
with the reference vector j^'j)'^"^. Based on the expan- 
sion 

= E ® I^bUb))}', (8) 



by summing over the reference subspace A for a fixe d 
value of js, one defines the reduced density matrix [36|: 

kA 

By construction, the density matrix is block-diagonal 
in both jB- In the warm-up phase, the reference sub- 
space becomes the 'medium' for the 'system' part in the 
B subspace. 

Truncation in the system sector is dictated by the 
density matrix. In standard DMRG applications for 
Hermitian problems where the eigenvalues of p are real 
non-negative, only the eigenvectors corresponding to the 
largest eigenvalues are kept during the DMRG process. 
Within the metric defining the Berggren ensemble, the 
GSM density matrix is complex-symmetric and its eigen- 
values are, in general, complex. In the straightforward 
generalization of the DMRG algorithm to the complex- 
symmetric case [2Sj| , one retains at most N^^^ eigenstates 
of/5^, 

P^{jB)\a)B ^'Wa\a)B, (10) 

having the largest nonzero values of \wa \ ■ Due to the nor- 
malization of l^j), the sum of all (complex) eigenvalues 
of /5b is equal to 1: 

Tr(p^) = ^^„ = l, (11) 

a 

i.e., the imaginary part of the trace vanishes exactly. 

Expressing the eigenstates in terms of the vectors 
\i)B hi B: 

i 

all matrix elements of the suboperators in these opti- 
mized states, 

B{a\0\P)B^Y.'^'^4 b{W)b, (13) 

are recalculated and stored. 

The warm-up procedure continues by adding to the 
system part the configurations containing particles in the 
second scattering shell {lj)2- As in the first step, one con- 
structs all many-body states {(0)2 } within this new shell 
and calculates corresponding matrix elements of subop- 
erators ([7]). The new vectors |i)b in the system sector 
are then obtained by coupling the states |a)s calculated 
in the first step with the vectors \{{lj)2})- 

Following the same prescription as before, one con- 
structs the set {kA in}'' of states coupled to in 
which the Hamiltonian is diagonalized. As in the pre- 
vious step, the new target state for the calculation of 
the reduced density matrix pB is defined as the one with 
maximum overlap with the reference state. Again, at 
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Warm-up Phase of GSM+DMRG 



{OCfi} 



FIG. 2: Schematic illustration of the GSM+DMRG proce- 
dure during the s*'' step of the warm-up phase. States {/ca} 
from A, previously optimized states ub, and states {{lj)s} 
constructed by occupying the s*'' shell with n particles are 
coupled to generate the new set of states {kA (E) is}'^ ~ 

{kA ^ {as ^ {ij):}y . 



most N^^pl eigenvectors of pB are retained and all matrix 
elements of suboperators for these optimized states are 
recalculated. This procedure, illustrated schematically in 
Fig. [21 continues until the last shell in B is reached, pro- 
viding a first guess for the wave function of the system 
in the whole ensemble of shells. At this point, all s.p. 
states have been considered, and all suboperators of the 
Hamiltonian H acting on states saved after truncation in 
B have been computed and stored. The warm-up phase 
ends and the so-called sweeping phase begins. 



B. Sweeping phase of GSM+DMRG 

Starting from the last scattering shell {lj)iast, the pro- 
cedure continues in the reverse direction (the 'sweep- 
down' phase) using the previously stored information. At 
this stage, the meaning of the medium and system parts 
changes as compared to the warm-up phase. 

In the sweeping phase, the states Ifc^) of the reference 
subspace A and the states \iprev) generated in the warm- 
up phase form the medium. The corresponding basis is: 



(14) 



The system part is generated by adding the scattering 
shells one at a time. 

The sweep-down process begins by constructing all 
possible states \i) from the shell {lj)iast and calculat- 
ing the corresponding suboperators of H . A new basis 
coupled to J'^ is then formed by coupling states 



with \i{jB))- The representation of H in this basis is con- 
Sweeping Phase of GSM+DMRG 
Sweep-down Phase Sweep-up Phase 



FIG. 3; Schematic illustration of the coupling between 
different configurations during the sweeping phase of 
GSM+DMRG. Left: sweep-down phase. Here, configurations 
{kA} axe coupled with previously generated configurations 
{iprev} and with {i} states. Right: sweep-up phase. Here, 
configurations {^a} are coupled with {i} and with configura- 
tions {iprev} in B generated in a previous sweep-down phase. 



structed using the Wigner-Eckart theorem by coupling 
suboperators acting in A, {iprev}, and {i} (the set of 
states As before, the target state 



I*,/) = 



E 



Jpart,jE 



{\^ Jp.rtik, iprev)) ® \iijB))}'\l5) 



is identified by picking up the eigenstate of H having 
the largest overlap with the reference state |\E'j)(°^ The 



density matrix is then constructed 



(16) 



and diagonalized for each value of jb- At this point, 
the truncation can be done in two different ways. In 
the first truncation method (i) at most iVopt eigenvec- 
tors of the density matrix with the largest nonzero mod- 
uli of eigenvalues are kept. This is precisely the trunca- 
tion technique that has been employed in the warm-up 
phase. The actual number of states retained may vary 
since one considers only eigenvectors with nonzero eigen- 
values. The second method (ii), based on the identity 
PT|) . is a gen eralization of the dynamical block selection 
approach [33| . Here we focus on controlling the numeri- 
cal error by selecting in each step of the procedure p, Np 
vectors with the largest moduli of the eigenvalues so that 
the condition 




< e 



(17) 



is satisfied. The quantity e in (jl7p can be viewed as 
the truncation error of the reduced density matrix. It 
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is worth noting that while the trace of the reduced den- 
sity matrix is strictly equal to one (jlip . this is no longer 
the case for the restricted sum of eigenvalues in (fT7|) . 
In particular, the real part of the reduced trace may be 
greater than one and the imaginary part may be nonzero. 
For that reason, in Eq. (fT7|) one considers the real part 
of the partial trace. The smaller e, the larger number 
Np of eigenvectors must be kept. In particular, for e=0, 
all eigenvectors with non-zero eigenvalues are retained. 
One should emphasize that Np may change from one step 
to another. Section IIVI discusses the convergence of the 
GSM-I-DMRG procedure with respect to A'opt and e. 

The matrix elements (|13l) in eigenvectors \a)B saved 
after the truncation are recalculated and stored. The 
procedure continues by adding the next shells one by 
one until the first scattering shell is reached. At each 
step during the sweep-down phase, all suboperators of H 
are stored. The sweep-down phase of GSM-f DMRG is 
schematically illustrated in the left portion of Fig. [31 

At this point, the procedure is reversed, and a sweep 
in the upward direction (the 'sweep-up' phase) begins. 
Using the information previously stored, a first shell is 
added, then a second one, etc. (see Fig. [3l right panel). 
The medium now consists of states in the reference sub- 
space A and states {iprev} in B that were generated 
during a previous sweep-down phase. The sweeping se- 
quences continue until convergence for target eigenvalue 
is achieved. 



IV. APPLICATIONS OF THE GSM+DMRG 
METHOD 



This work describes the first GSM-I-DMRG treatment 
of open-shell proton-neutron nuclei. As illustrative ex- 
amples, we take ^Li and *Li described schematically as 
interacting nucleons outside the closed core of '*He. The 
neutron one-body potential in Eq. ([2]) is a Woods-Saxon 
(WS) potential with radius i?o = 2 fm and diffuseness 
d = 0.65 fm. The spin-orbit strength V^o — 7.5 MeV and 
the depth of the central potential Vb=47 MeV are fixed to 
reproduce the experimental energies and widths of 3/2^ 
and 1/2 j" resonances in ^He. For the protons, the same 
WS average potential is supplemented by the Coulomb 
potential generated by a uniformly charged sphere of ra- 
dius i?o and charge Q—+2e. 

The two-body interaction in ([2]) is represented by a 
finite-range Surface Gaussian Interaction (SGI) [s^: 



neutrons Nyi 



J,T 



= Va{J,T)exp 



X <5(|ri| + |r2|-2i?o). 



(18) 



Vr 



3,0 



"30 



-1), 

/33o(A^. -1), 



(19) 



The strengths Vq{J,T) are the same as in Ref. [33 ■ The 
T=0 couplings depend linearly on the number of valence 



where aio = -600 MeV ■ fm^ (3iq = -50 MeV ■ fm^ 
ago = -625 MeV • fm^ and ^30 = -100 MeV • fm^ 

The above set of parameters has not been optimized 
to reproduce the actual structure of ^Li and *Li. Our 
choice of interaction is motivated by the fact that the 
main purpose of this study is to test the DMRG pro- 
cedure for proton-neutron systems for which the exact 
GSM diagonalization is still possible. In this context, 
"our" ^Li and ®Li systems should be viewed as three- 
and four-nucleon cases, respectively. 

Following the method described in Ref. [Ill, s.p. bases 
for protons and neutrons are generated by their respec- 
tive spherical Hartree-Fock (HF) potentials correspond- 
ing to the GSM Hamiltonian ([2]). Neutron and proton 
valence spaces include the 0^3/2 HF poles as well as the 
scattering shells {p3/2}c in the complex fc-plane, and the 
{Pi/2}c, si/2, and ^5/2 real-energy continua. 

The Lt.„ (Lt-„ ) contour, along which the scat- 
tering 1^^3/2 (7rp3/2) shells are distributed, is defined 
by a triangle with vertices at: (i?e(fc), /m(/c)) = 
(0.0, 0.0), (0.17, -0.1), (2, 0.0), and a segment along the 
i?e(fc)-axis from (2.0,0.0) to (8.0,0.0) (in units of fm'^). 
Each segment of these contours is discretized with two 
points corresponding to the abscissa of the Gauss- 
Legendre quadrature. Hence, we take 6 non-resonant 
continuum shells from {^3/2)0 both for protons and neu- 
trons. The real-energy pi/2 continuum shells are dis- 
tributed along the segment [(0,0), (8,0)] which is dis- 
cretized with 6 points for protons and neutrons. The 
real-energy Si/2 and d^^2 continua arc included in the 
valence space as well. They are distributed along the 
real /c-axis along the segment [(0,0), (3,0)]. We take six 
{si/2}c and {dQ/2}c discretization points, both for pro- 
tons and neutrons. The Osi/2 poles are not included in 
the valence space as they are assumed to be occupied in 
the core of ^He. The total number of shells A^sh in the 
GSM configuration space is then equal to 50. 



The three- nucleon case: =3/2 ground state 
of ^Li 



The s.p. basis of ^Li is generated by the HF potential 
(calculated separately for protons and neutrons and for 
each partial wave). It contains bound s.p. P3/2 states at 
energies —5.605 MeV (neutrons) and —7.098 MeV (pro- 
tons). These s.p. states generate the pole space in the 
many-body GSM framework and the reference subspace 
A in DMRG. As discussed above, the total number of res- 
onant (P3/2) and non-resonant ({p3/2}c, {pi/2}c, {si/2}c, 
{^5/2}c) shells for protons and neutrons is 50. The di- 
mension of the Lanczos space spanned by one valence 
proton and two valence neutrons in these 50 shells, i.e.. 
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the dimension of the GSM matrix, is D=7,796. The 
ground state energy E^^ = (-26.6620, 0.2486) MeV has 
a non- vanishing, unphysical imaginary part. This spu- 
rious width comes from the fact that the discretization 
along the contours is not precise enough to effectively 
fulfill the completeness relation ([T]). This problem will 
be addressed in Sec. IIV A3I by increasing the number of 
points along the contour. In what follows, we shall study 
the convergence of the GSM+DMRG method by varying 
either the number of eigenvectors A^opt kept during the 
sweeping phase (see Sec. IIV A ip . or the precision e of 
the density matrix (see Sec. IIV A2p . 
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1. DMRG truncation with fixed A^opt 

The number of eigenvectors of pB with the largest 
nonzero moduli of eigenvalues kept at each iteration dur- 
ing the warm-up phase is limited to nI^^.—2Q. This num- 
ber corresponds to the total number of states {n;jB} in 
the subspace B that can be coupled with states in A to 
yield configurations with J'^ = 3/2~. 

The actual number of eigenvectors kept in the warm- 
up phase may be less than N^^^ since most of eigenvec- 
tors have vanishing eigenvalues. The non-resonant con- 
tinuum shells involved in B are ordered according to the 
sequence: 



FIG. 4: The ground-state energy of ■'Li as a function of the 
step number in the sweeping phase, Nstep- For both values of 
A'^opt, the maximum number of eigenvectors kept during the 
warm-up phase is Ai'^pt=26. The step zero corresponds to the 
result obtained at the end of the warm-up phase. The solid 
line marks the exact GSM value obtained by means of the 
direct Lanczos diagonalization. 
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where index i denotes the position of scattering shells on 
their respective contours, beginning with those closest to 
the A;=0 origin. 



Figure [4] illustrates the convergence of the 
GSM-I-DMRG procedure with respect to the step 
number in the sweeping phase for A'opt=40 and 60. The 
results are identical for both values of A"opt, and the 
DMRG ground-state energy converges to the value of 
i?e(£;DMRG) = -15.176 MeV. This exceeds by -7 MeV 
the exact GSM value of i?e(i;cx)=-22.662MeV, ob- 
tained by the direct Lanczos diagonalization of the GSM 
Hamiltonian. 

Clearly, when applied to the case shown in Fig. 01 the 
GSM-(-DMRG procedure breaks down. The reason for 
this failure is not related to a too small value of Nopt- in- 
deed, in the case A'opt = 60 the largest number of eigen- 
vectors of the density matrix with non-zero eigenvalues 
is equal to 50. The GSM-hDMRG iterative method is 
trapped in a local minimum, a not uncommon feature 
of the standard DMRG procedure. A further increase 
of A^opt does not change the final results which is fully 
converged. To understand the origin of the failure, let us 
analyze the GSM-f DMRG wave function in some detail. 



TABLE I: Real part of the squared shell-model amplitudes in 
the J'^ = 3/2~ ground state wave function of '^Li obtained 
in GSM-I-DMRG (A'^opt = 60) and through the exact Lanczos 
diagonalization . 

To this end, the J'^ = 3/2^ g.s. wave function of ^Li is 
decomposed as follows: 

I*) = Cp3\p^) + Cs2p\s'^p) + Cd2p\(fp) + Cspd\spd), (21) 

where c^s are the amplitudes associated with different 
three-nucleon GSM configurations The (real parts) 
of squared amplitudes are shown in Table U for the 
GSM-f GDMRG wave function corresponding to Aopt=60 
and for the exact GSM wave function. As compared to 
the exact result, the \p^) parentage amplitude is overes- 
timated and the \spd) component is totally absent in the 
GSM-I-DMRG wave function. The latter can be under- 
stood by observing that (i) only shells with 1=1 span the 
reference subspace and (ii) during the GSM-I-DMRG 
procedure, scattering shells are added one by one. Conse- 
quently, when the first positive-parity shell (in our case, a 
■KSi/2 shell, see (PU]) ) is added, the \spd) component can- 
not be generated as the first c?5/2 shell is added only later. 
When the first non-resonant shell is included, the 
I spd) configuration cannot be generated either, because 
states with one particle in previously considered s-shells 
are not kept in the process of optimization due to the par- 
ity conservation. Therefore, the \spd) configuration never 
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enters the DMRG wave function; hence, GSM+DMRG 
converges to a wrong solution. To prevent this patho- 





FIG. 5: Schematic illustration of the reference subspace A' 
during the warm-up phase of GSM+DMRG for J"" = 3/2" 
wave functions of ^Li. In addition to the 7i"(0p3/2), 1^(0^3/2), 
7r(0pi/2)j and 1^(0^1/2) poles, two scattering shells 1^(1^/2 and 
are now included to generate \spd) configurations during 
the DMRG procedure. See text for more details. 



logical behavior, we add to the reference subspace A two 
positive-parity scattering shells 7rd5/2 and to form 

a new reference subspace A' (see Fig. [5]). We arbitrarily 
choose the last T^d^j^ and shells in the sequence 

(1^0)1 . (As we shall see later, any other positive parity 
shells can be chosen as well.) The role played by the ad- 
ditional positive parity shells is to generate missing SM 
couplings in the wave function. The new reference sub- 
space A' is used for the construction of the set {fc^/ (8) is} 
and the density matrix . At each iteration during the 
warm-up phase, the density matrix contains the corre- 
lations due to the additional positive-parity orbits. In 
this way we assure that no possible couplings are missing 
during the warm-up phase. We use the same reference 
state as before, i.e., (which is only generated by the 
resonant shells) to select the target state among the 
eigenstates of H. By the time the first sweep starts, the 
two shells (7r(i5/2, 1^^5/2) are included in B and the pro- 
cedure is carried out as was described in Sec. IIIII with the 
reference subspace A. 



Real and imaginary parts of the ground state energy, 
obtained using the extended reference subspace A' of 
Fig. [5l are plotted in Fig. [6] for different values of A"opt- 
For A'opt=40, one can see pronounced quasi-periodic os- 
cillations in both real and imaginary parts of the energy. 
These oscillations have the periodicity of 96 steps cor- 
responding to two consecutive sweeps: sweep-down and 
sweep-up, each consisting of 48 steps. The energy oscilla- 
tions rapidly diminish with increasing N^pt, and the cal- 
culated energy E'dmrg converges to the GSM benchmark 
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FIG. 6: The GSM-I-DMRG ground state energy (real part, 
top; imaginary part, bottom) of ^Li for Afopt=40, 60, and 80 as 
a function of the step number during the sweeping phase. Two 
non-resonant continuum shells (ti" 6(5/2, ^^^5/2) are included in 



the reference subspace during the warm-up phase {N'^^^^—2&) 



resuh: E^^ = (-26.6620, 0.2486) MeV. For Afopt=80, the 
deviation from the benchmark result is less than 1 keV 
for the real part of the energy, and less than 0.1 keV for 
the corresponding imaginary part. 

The rank d^'^^ of largest matrix to be diagonalized 
grows almost linearly with iVopt, from c?^^''''=716 for 
TV^pt = 40 (- 9.1 % of the dimension D of the GSM 
matrix) to '^'^=1469 (~ 19 % of D) for iVopt=80. One 
should keep in mind that d^'^^ is almost independent 
of the continuum discretization density [HI, i.e., the 
number of scattering shells considered. Hence, the ra- 



tio dg^'^/D decreases rapidly with the number of valence 
shells 29]. 

The GSM-I-DMRG energy averaged over steps of the 
fourth sweep, as well as the minimum and maximum 
energy value reached during this sweep, are plotted in 
Fig. [7] for various iVopt. As A^opt increases, the am- 
plitude of energy, defined as the difference between the 
maximum and the minimum of real (top) or imaginary 
(bottom) part of GSM-f DMRG energy during the fourth 
sweep, decreases monotonously. Moreover, with increas- 
ing A"opt, both real and imaginary parts of the average 
energy converge exponentially to the exact value. Re- 
sults of x^-analysis are shown by solid lines in Fig. [T] 
Asymptotic values extracted in this way, Rg{E^^^^q) — 
-26.6622 ± 0.0002 MeV and Im{E'^^f^) = 0.248580 ± 
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0.246 



FIG. 7: Average value of the ground state energy (top: real 
part; bottom: imaginary part) of ^Li in GSM+DMRG for 
different values of A^'opt- Error bars correspond to the dif- 
ference between extremum values of -Re(-E'DMRG) (top) and 
/m(i5DMRG) (bottom) at the fourth sweep. Two scattering 
shells {jrd^/2,vd5/2) are included in the reference subspace 
during the warm-up phase {nI^^—2&). 



Nopf 


|-£'max 1 


fle(S„,ax) 


Jm(Smax) 


i?e(Savo) 


/m(i5avc) 


40 


22.6489 


-22.6475 


0.2470 


-22.5270 


0.2468 


50 


22.6605 


-22.6591 


0.2484 


-22.61844 


0.2478 


60 


22.6631 


-22.6617 


0.2485 


-22.6510 


0.2484 


70 


22.6634 


-22.6620 


0.2486 


22.6609 


0.2486 


80 


22.6634 


-22.6620 


0.2486 


-22.6619 


0.2486 



TABLE II: Modulus, real and imaginary part of -Bmax de- 
fined as the DMRG energy with the greatest modulus during 
the last sweep. The real and imaginary parts of the average 
energy -Bavc at the fourth sweep (corresponding to the case 
presented in Fig [7] are also shown for comparison. 
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FIG. 8: Real part of the ground state energy of ^Li as a 
function of Aopt. For each Aopt, the energy with the largest 
modulus during the last sweep iSmax is selected. Two scatter- 
ing shells {nd^p, 2^^5/2) are included in the reference subspace 
during the warm-up phase (A^pj=26). 



0.000004 MeV, reproduce the exact GSM result very well. 
The feature of an exponential convergence of the step- 
averaged GSM-I-DMRG energies may be useful when es- 
timating eigenvalues based on results obtained with rel- 
atively small iVopt. 

To illustrate how the generalized variational principle 
^ works, let us consider the energy with the great- 
est modulus, i?max, Calculated in DMRG during the last 
sweep. The values of E„iax and i?avc, the energy averaged 
during the last sweep (corresponding to Fig[7]), are shown 
in Table mi for different iVopt. One can clearly see that the 
closer the wave function calculated with DMRG is to the 
exact wave function as iVopt increases, the larger |i?max| 
is. Hence, in this case, the modulus of energy reaches a 
maximum at the local extremum of the functional E[^] 
corresponding to the ground state energy of ''Li. 

The convergence to the exact value is faster by con- 
sidering iSmax for each truncation A'opt instead of se- 
lecting the average value i?avc (cf. Table [n| . The real 
part i?e(£'niax) converges exponentially to R^{E^^^f^) = 
-26.6621±0.0002MeV while the imaginary part of E^^^ 
doesn't follow the exponential behavior. 

The choice of positive-parity scattering shells to be in- 
cluded in the reference subspace is somehow arbitrary. 



The only important point is that by including both posi- 
tive and negative parity shells during the warm-up phase, 
one can generate many-body configurations that would 
not appear otherwise. For that reason, one can replace 
d5/2 with si/2 scattering shells without changing the out- 
come of the GSM-I-DMRG procedure. To illustrate this. 
Fig. [Ushows the GSM-I-DMRG results with the extended 
reference subspace A' containing either two (7rsi/2, 1^51/2) 
or scattering shells. It is seen that the con- 

verged value of the GSM-f DMRG energy is the same in 
both cases. 

A different way to generate the missing components of 
the wave function is to demand that at each step during 
the warm-up phase at least one state from each {n;jB} 
family is kept after truncation. We take up to n'^^ eigen- 
vectors of the density matrix with largest nonzero eigen- 
values, where A^^^pt is equal to the number of different 
families js} which contribute to the GSM-I-DMRG 
wave function. If certain families are not represented 
in this set of eigenvectors, we add one state for each 
such family even if the corresponding eigenvalue equals 
zero. Hence, the actual number of vectors kept during 
the warm-up phase almost always exceeds A^opt- Using 
this additional condition, one may employ a standard 



9 



T 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 f 




I I I I I I I I I I I I I I I I I I L_ 

20 40 60 80 100 120 140 160 180 

Nstep 



FIG. 9: The real part of the ground state energy of "^Li as a 
function of A'stcp. The sohd hue shows resuhs obtained using 
an extended reference space spanned on the {p3/2,Pi/2) poles 
space and two (7rsi/2, i^Si/2) scattering shells. Results de- 
picted by squares have been obtained using a standard setup 
for an extended reference subspace as in Fig. |6l Open circles 
show results obtained with a reference subspace spanned on 
the (p3/2,Pi/2) poles and demanding that at least one state 
from each {n; js} family is kept during the warm-up phase. 
The values of N^ll^26 and No pt=80 were used in all cases. 
See text for more details. 



setup for the reference subspace (i.e., A is spanned by s.p. 
poles) . Results using this GSM-I-DMRG strategy are also 
shown in Fig. [9l The minimal number of states which 
are kept in the warm-up phase is, in this case, Nj^^^=26. 
In spite of a rather different energy at the beginning of 
the sweeping phase, the exact GSM-I-DMRG energy is 
reproduced. Moreover, the use of an extended reference 
subspace improves convergence. The rank of the largest 
matrix to be diagonalized in this case, dg'''^=1469, is in- 
dependent of the algorithm chosen. 



2. Truncation governed by the trace of the reduced density 
matrix 

In the examples described in Sec. lIV A !( the maximum 
number of states iVopt is kept fixed at each step of the 
sweeping phase. This does not mean that the number 
of eigenvectors retained in the sweeping phase is always 
constant or equal iVopt; only the eigenvectors of the den- 
sity matrix with no n- vanishing eigenvalues are kept. In 
this section, we shall investigate the GSM-I-DMRG algo- 
rithm in which the number of states Np kept at any step 
in the sweeping phase depends on the condition (jl7|) for 
the trace of the density matrix. The real part of the 
J'^ — 3/2 j~ eigenvalue in ^Li is shown in Fig. [10] for sev- 
eral values of e. As in the previous examples, the ref- 
erence subspace A' is spanned by the HF poles and two 
scattering shells (7rc?5/2, 1^^5/2)- In the warm-up phase, 

we keep iV^pt eigenvectors of the density matrix (up to 
26) and additionally require that at least one state of each 



{n; Jb} family is retained. As before, N^p{ is equal to the 
number of different {?t.;jb} families so the total number 
of eigenvectors kept at each iteration step is greater or 
equal to iV^p^. For low-precision calculations (e = 10"'*), 
the resulting energy oscillates and approaches a value 
which deviates from the correct result by ~1.9MeV. The 
amplitude of oscillations as a function of A'step quickly 
decreases with decreasing e. For e = 10~^, the precision 
of the converged GSM-I-DMRG energy value is ~ 0.2 keV 
for both real and imaginary parts. 
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FIG. 10; The ground state energy (real part) of ^Li as a 
function of the step number during the sweeping phase for 
five values of the truncation error e of the reduced density 
matrix, see Eq. (I17p . Two scattering shells (716^5/2,1^^5/2) are 
included in the reference subspace during the warm-up phase. 

Obviously, the dimension of the largest matrix to be 
diagonalized depends on the required truncation error e. 
In the studied case, d^^^ changes from 273 for e = 10"''' 
to 1327 for e = 10~* with the average number of vectors 
kept during the sweeping phase increasing from ~15 to 
~46. In the truncation scenario with fixed A'opt, the 
number of saved vectors, averaged over one sweep, is 
^--^59 for iVopt=80. In general, for the same precision of 
GSM-I-DMRG energies, the average number of vectors 
kept during the sweeping phase is smaller if the trunca- 
tion is done dynamically according to the trace of density 
matrix than by fixing the maximum number of eigenvec- 
tors iVopt. 

The GSM-I-DMRG energy averaged over steps of the 
third sweep as well as the minimum and maximum energy 
reached during this sweep are plotted in Fig. [TT] as a 
function of e. The GSM+DMRG error, i.e., the energy 
difference with respect to the exact GSM result, decreases 
fast with decreasing e. The real and imaginary parts of 
DMRG energy satisfy to a good approximation the power 
law 

Re (Eumrg) - Re (E,^) + ae^ (22) 
proposed in Ref. [s^ to control the accuracy of the 
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FIG. 11: The ground state average energy (top: real part; 
bottom: imaginary part) of ^Li as a function of the trunca- 
tion error e of the reduced density matrix. Two scattering 
shells (7rd5/2, 1^^5/2) are included in the reference subspace 
during the warm-up phase. Error bars correspond to a dif- 
ference between the extremum values of Re{EDMRG) at the 
third sweep. The solid line shows the results of analysis, 
assuming power-law convergence (|22p . 



e. 


|-£'max 1 








/m(i5avc) 


10" 


4 


20.9221 


-20.9209 


0.2240 


-20.7751 


0.2157 


10^^ 


5 


22.5575 


-22.5562 


0.2416 


-22.4870 


0.2434 


10" 


6 


22.6532 


-22.6519 


0.2485 


-22.6474 


0.2479 


10" 


7 


22.6621 


-22.6607 


0.2481 


-22.6602 


0.2483 


10" 


8 


22.6632 


-22.6618 


0.2486 


-22.6618 


0.2484 



TABLE III: Modulus, real and imaginary part of -Emax defined 
as the energy having the greatest modulus during the last 
sweep. The real and imaginary parts of the average energy 
Eaverage at the fourth sweep (see Fig fTTj) are also shown for 
comparison. 
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FIG. 12: Real part of the ground state energy of ^Li as a 
function of e. At each e, the energy is selected according to 
the greatest modulus during the last sweep. 



DMRG method for Hermitian problems. The results of 
a x^"fit to Eq. ([221) are shown in Fig. [TT] The asymp- 
totic values extracted in this way are ReiE^^^f^) — 
-26.66192±2 ■ lO''^ MeV and Im{El^^Q) = 0.24844 ± 
3 • 10~^ MeV and agree very well with the exact GSM 
energy. 

In Table IIIII we compare the average complex energy 
i?avo at the third sweep for different values of e (cor- 
responding to Fig [11] and the complex energy -Emax 
(the energy with the greatest modulus during the last 
sweep). As in the previous case where Nopt was fixed, 
the modulus of -Emax reaches a maximum when E-a^ay, 
is equal to the exact GSM energy. As can be seen in 
Fig. [T21 the real part of -Emax exhibits a power-law be- 
havior with an extrapolated value equal to Rg{E^^^^q) = 
-26.660.5 ±0.0012 MeV. 



3. Treatment of spurious width 

As mentioned in Sec. IIV Al the imaginary part of 
the GSM energy Eox = (-26.6620, 0.2486) MeV is non- 



physical for it has a negative width. This is due to the 
fact that the contour discretization is not sufficiently pre- 
cise to guarantee the completeness relation ([ij . This spu- 
riosity can be taken care of by increasing the number of 
points along the integration contour. In the largest calcu- 
lation we have done for the ground state of ''Li, we took 
67 points along the contour L+p^^^, 24 along L+p^^^, and 

12 points along the contours Lt.^ and L't j . The 
neutron valence space is the same as the proton space 
except for the contour L'^.^^^^ where 74 points are consid- 
ered. The model space corresponds to 239 shells and the 
dimension of the ground state = 3/2^ GSM Hamil- 
tonian matrix is 1,459,728. 

In order to perform calculations within this huge va- 
lence space, we have developed a parallel version of the 
DMRG code. At each step during the DMRG procedure, 
calculations of the matrix elements of the suboperators 
([7]) and hamiltonian ^ are distributed among the pro- 
cessors. Our calculations were carried out on the CRAY 
XT4 Jaguar supercomputer at the Oak Ridge National 
Laboratory. 

The real part of the ground state energy and the fit 
according to the relation ([22]) are plotted in Fig [131 For 
each e the energy Emax with the greatest modulus during 
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the third sweep is considered. The extrapolated value is 
■SeC^^DMRG) = -21.6834 ± O.OOlOMeV. The real part of 
Emax for e = 5-10-1" is -21.6820 MeV and the amplitude 
during the last sweep is 2.275 keV; hence, convergence 
has almost been reached. Here, the largest matrix has 
a dimension 3,348. The imaginary part (which does not 
follow the power law behavior) varies from 0.00100 MeV 
at e = lO-'^ to 0.00075 MeV at e = 5 • 10"^" (its ampli- 
tude during the last sweep is 0.065 keV). This example 
nicely demonstrates the validity of the many-body com- 
pleteness relation in GSM. 
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FIG. 13: Real part of the ground state energy of ^Li in a 
model space made of 239 shells. At each e, the energy is 
selected according to the greatest modulus during the third 
sweep. See text for details. 



set of eigenvectors, we add one state for each such fam- 
ily even if the corresponding eigenvalue equals zero. In 
the sweeping phase, we follow the truncation strategy of 
Sec. lIVA2l 




N 



B. The four-nucleon case: J'" = 2^ ground state of 

The GSM Hamiltonian for Li is the same as for ^Li, 
the only modification being the change of T=0 cou- 
plings in Eq. (jl9p due to the different number of neu- 
trons Nj;=3. The HF procedure yields two bound s.p. 
states: Cp.^^^ = -8.556 MeV and Cp^^^ = -12.788 MeV, 
for neutrons and protons, respectively. Shells of the non- 
resonant continuum are distributed in the complex k- 
plane using the same contours and the same discretiza- 
tion scheme as in the ^Li case. The dimension of the 
Lanczos space spanned by one valence proton and three 
valence neutrons in 50 shells, i.e., the dimension of the 
GSM matrix in ^Li, is D=170,198. The GSM+DMRG 
results presented in this section are obtained using the 
truncation criterion ([T7]). As discussed in Sec. IIV A 21 
this criterion is somewhat more efficient than the condi- 
tion based on fixing the maximum number of eigenvectors 

The truncation method employed in the warm-up 
phase follows that of Sec. IIV A II The reference sub- 
space A' is spanned by the pole states and two scatter- 
ing shells {Trd^/2,^d^/2)- We take up to 50 eigenvectors 
of the density matrix with the largest nonzero eigenval- 
ues. If certain {n; Jb} families are not represented in this 



FIG. 14: The real part of the ground state energy of Li as a 
function of step number during the sweeping phase. Results 
are shown for different values of the truncation error e. Two 
non-resonant continuum shells {Tid^i2,vd^/2) are included in 
the reference subspace during the warm-up phase. 

Figure [H shows the DMRG-I-GSM results for ^Li for 
three values of e. The exact energy of the 2^ state ob- 
tained by the direct Lanczos diagonalization of the GSM 
matrix is E'ox = (-19.19451, 0.I336I) MeV. The corre- 
sponding energies averaged over the third sweep are plot- 
ted in Fig. [151 For e = 10"^, the largest matrix to be di- 
agonalized has a rank ^^==1,446 (~ 0.8% of For e = 
10-^ one obtains Ey^mko = (-19.19415, 0.13355)MeV, 
i.e. the real part of the GSM-I-DMRG energy devi- 
ates only by 0.4 keV from the exact value, while for the 
imaginary part, the deviation is less than 0.06 keV. The 
largest matrix to be diagonalizcd in this case has a rank 
dg^''=20,535 (~ 12 % of D). 

V. CONCLUSIONS 

This work describes the first application of the DMRG 
method to two-fluid, open many-fermion systems rep- 
resented by complex-symmetric Hamiltonians. Calcu- 
lations were carried out for proton-neutron systems ^Li 
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FIG. 15: Similar as in Fig. [11] except for the ground state 
average energy of *Li. 



(three-nucleon problem) and ^Li (four-nucleon problem). 
As compared to our previous work [29|, two significant 
improvements of the GSM+DMRG technique have been 
made. The first improvement concerns the recognition of 
the appropriate target state in the warm-up phase. The 
second development relates to the truncation strategy in 
the sweeping phase. 

There are situations in which the DMRG procedure 
yields a fully converged but incorrect solution. In or- 
der to understand and prevent this pathological behav- 
ior, we studied a necessary and sufficient condition for 
the GSM-f DMRG method to yield a correct eigenvalue. 
The essential condition is to assure that all possible cou- 
plings in the many-body wave function, allowed by the 
symmetries of the problem and the configuration space, 
are present in the warm-up phase. We propose different 
strategies to guarantee this crucial requirement. 

Two truncation schemes for the number of retained 



vectors in the sweeping phase of DMRG were investi- 
gated: the fixed-A'^opt scheme (Sec. IIV A 1|) and the dy- 
namic truncation (Sec. IIV A 2^ . We conclude that the 
two strategies are to a large extent equivalent; they both 
exhibit the excellent convergence properties to the bench- 
mark GSM result. In both cases, one finds the quasi- 
periodic oscillations of GSM-I-DMRG energy as a func- 
tion of A'stcp with extensive plateaux. 

The GSM-f DMRG energy averaged over one sweep ex- 
hibits excellent exponential convergence with A'opt which 
allows to deduce the asymptotic value with good preci- 
sion. Also, i?DMRG exhibits excellent convergence as a 
function of the truncation error e. This feature makes it 
possible to control the accuracy of GSM+DMRG calcu- 
lations. 

The dynamic truncation, fixing a condition on the 
trace of the reduced density matrix, yields results of sim- 
ilar accuracy for ''Li and *Li, i.e., for systems having very 
different configuration spaces. This offers a possibility to 
compare the convergence in different quantal systems at 
the same value of e. 

The encouraging features of the proposed 
GSM-I-DMRG approach open the window for sys- 
tematic and high-precision studies of complex, weakly 
bound nuclei, such as halo systems, which require large 
configuration spaces involving s.p. states of different 
parities (both in the pole space and in the scattering 
space). Generally, the improvements of the DMRG 
approach proposed in this work can be of interest in the 
context of other multiparticle open quantum systems, 
as well as for other DMRG calculations involving 
non-Hermitian Hamiltonians. 
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